home *** CD-ROM | disk | FTP | other *** search
/ Game.EXE 2001 January / Game.EXE_01_2001.iso / demos / Blade of Darkness / data1.cab / Program_Executable_Files / Lib / PythonLib / random.py < prev    next >
Encoding:
Python Source  |  2000-11-16  |  7.2 KB  |  287 lines
  1. #    R A N D O M   V A R I A B L E   G E N E R A T O R S
  2. #
  3. #    distributions on the real line:
  4. #    ------------------------------
  5. #           normal (Gaussian)
  6. #           lognormal
  7. #           negative exponential
  8. #           gamma
  9. #           beta
  10. #
  11. #    distributions on the circle (angles 0 to 2pi)
  12. #    ---------------------------------------------
  13. #           circular uniform
  14. #           von Mises
  15.  
  16. # Translated from anonymously contributed C/C++ source.
  17.  
  18. from whrandom import random, uniform, randint, choice # Also for export!
  19. from math import log, exp, pi, e, sqrt, acos, cos, sin
  20.  
  21. # Housekeeping function to verify that magic constants have been
  22. # computed correctly
  23.  
  24. def verify(name, expected):
  25.     computed = eval(name)
  26.     if abs(computed - expected) > 1e-7:
  27.         raise ValueError, \
  28.   'computed value for %s deviates too much (computed %g, expected %g)' % \
  29.   (name, computed, expected)
  30.  
  31. # -------------------- normal distribution --------------------
  32.  
  33. NV_MAGICCONST = 4*exp(-0.5)/sqrt(2.0)
  34. verify('NV_MAGICCONST', 1.71552776992141)
  35. def normalvariate(mu, sigma):
  36.     # mu = mean, sigma = standard deviation
  37.  
  38.     # Uses Kinderman and Monahan method. Reference: Kinderman,
  39.     # A.J. and Monahan, J.F., "Computer generation of random
  40.     # variables using the ratio of uniform deviates", ACM Trans
  41.     # Math Software, 3, (1977), pp257-260.
  42.  
  43.     while 1:
  44.         u1 = random()
  45.         u2 = random()
  46.         z = NV_MAGICCONST*(u1-0.5)/u2
  47.         zz = z*z/4.0
  48.         if zz <= -log(u2):
  49.             break
  50.     return mu+z*sigma
  51.  
  52. # -------------------- lognormal distribution --------------------
  53.  
  54. def lognormvariate(mu, sigma):
  55.     return exp(normalvariate(mu, sigma))
  56.  
  57. # -------------------- circular uniform --------------------
  58.  
  59. def cunifvariate(mean, arc):
  60.     # mean: mean angle (in radians between 0 and pi)
  61.     # arc:  range of distribution (in radians between 0 and pi)
  62.  
  63.     return (mean + arc * (random() - 0.5)) % pi
  64.  
  65. # -------------------- exponential distribution --------------------
  66.  
  67. def expovariate(lambd):
  68.     # lambd: rate lambd = 1/mean
  69.     # ('lambda' is a Python reserved word)
  70.  
  71.     u = random()
  72.     while u <= 1e-7:
  73.         u = random()
  74.     return -log(u)/lambd
  75.  
  76. # -------------------- von Mises distribution --------------------
  77.  
  78. TWOPI = 2.0*pi
  79. verify('TWOPI', 6.28318530718)
  80.  
  81. def vonmisesvariate(mu, kappa):
  82.     # mu:    mean angle (in radians between 0 and 2*pi)
  83.     # kappa: concentration parameter kappa (>= 0)
  84.     # if kappa = 0 generate uniform random angle
  85.  
  86.     # Based upon an algorithm published in: Fisher, N.I.,
  87.     # "Statistical Analysis of Circular Data", Cambridge
  88.     # University Press, 1993.
  89.  
  90.     # Thanks to Magnus Kessler for a correction to the
  91.     # implementation of step 4.
  92.  
  93.     if kappa <= 1e-6:
  94.         return TWOPI * random()
  95.  
  96.     a = 1.0 + sqrt(1.0 + 4.0 * kappa * kappa)
  97.     b = (a - sqrt(2.0 * a))/(2.0 * kappa)
  98.     r = (1.0 + b * b)/(2.0 * b)
  99.  
  100.     while 1:
  101.         u1 = random()
  102.  
  103.         z = cos(pi * u1)
  104.         f = (1.0 + r * z)/(r + z)
  105.         c = kappa * (r - f)
  106.  
  107.         u2 = random()
  108.  
  109.         if not (u2 >= c * (2.0 - c) and u2 > c * exp(1.0 - c)):
  110.             break
  111.  
  112.     u3 = random()
  113.     if u3 > 0.5:
  114.         theta = (mu % TWOPI) + acos(f)
  115.     else:
  116.         theta = (mu % TWOPI) - acos(f)
  117.  
  118.     return theta
  119.  
  120. # -------------------- gamma distribution --------------------
  121.  
  122. LOG4 = log(4.0)
  123. verify('LOG4', 1.38629436111989)
  124.  
  125. def gammavariate(alpha, beta):
  126.         # beta times standard gamma
  127.     ainv = sqrt(2.0 * alpha - 1.0)
  128.     return beta * stdgamma(alpha, ainv, alpha - LOG4, alpha + ainv)
  129.  
  130. SG_MAGICCONST = 1.0 + log(4.5)
  131. verify('SG_MAGICCONST', 2.50407739677627)
  132.  
  133. def stdgamma(alpha, ainv, bbb, ccc):
  134.     # ainv = sqrt(2 * alpha - 1)
  135.     # bbb = alpha - log(4)
  136.     # ccc = alpha + ainv
  137.  
  138.     if alpha <= 0.0:
  139.         raise ValueError, 'stdgamma: alpha must be > 0.0'
  140.  
  141.     if alpha > 1.0:
  142.  
  143.         # Uses R.C.H. Cheng, "The generation of Gamma
  144.         # variables with non-integral shape parameters",
  145.         # Applied Statistics, (1977), 26, No. 1, p71-74
  146.  
  147.         while 1:
  148.             u1 = random()
  149.             u2 = random()
  150.             v = log(u1/(1.0-u1))/ainv
  151.             x = alpha*exp(v)
  152.             z = u1*u1*u2
  153.             r = bbb+ccc*v-x
  154.             if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= log(z):
  155.                 return x
  156.  
  157.     elif alpha == 1.0:
  158.         # expovariate(1)
  159.         u = random()
  160.         while u <= 1e-7:
  161.             u = random()
  162.         return -log(u)
  163.  
  164.     else:    # alpha is between 0 and 1 (exclusive)
  165.  
  166.         # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
  167.  
  168.         while 1:
  169.             u = random()
  170.             b = (e + alpha)/e
  171.             p = b*u
  172.             if p <= 1.0:
  173.                 x = pow(p, 1.0/alpha)
  174.             else:
  175.                 # p > 1
  176.                 x = -log((b-p)/alpha)
  177.             u1 = random()
  178.             if not (((p <= 1.0) and (u1 > exp(-x))) or
  179.                   ((p > 1)  and  (u1 > pow(x, alpha - 1.0)))):
  180.                 break
  181.         return x
  182.  
  183.  
  184. # -------------------- Gauss (faster alternative) --------------------
  185.  
  186. gauss_next = None
  187. def gauss(mu, sigma):
  188.  
  189.     # When x and y are two variables from [0, 1), uniformly
  190.     # distributed, then
  191.     #
  192.     #    cos(2*pi*x)*sqrt(-2*log(1-y))
  193.     #    sin(2*pi*x)*sqrt(-2*log(1-y))
  194.     #
  195.     # are two *independent* variables with normal distribution
  196.     # (mu = 0, sigma = 1).
  197.     # (Lambert Meertens)
  198.     # (corrected version; bug discovered by Mike Miller, fixed by LM)
  199.  
  200.     global gauss_next
  201.  
  202.     if gauss_next != None:
  203.         z = gauss_next
  204.         gauss_next = None
  205.     else:
  206.         x2pi = random() * TWOPI
  207.         g2rad = sqrt(-2.0 * log(1.0 - random()))
  208.         z = cos(x2pi) * g2rad
  209.         gauss_next = sin(x2pi) * g2rad
  210.  
  211.     return mu + z*sigma
  212.  
  213. # -------------------- beta --------------------
  214.  
  215. def betavariate(alpha, beta):
  216.  
  217.     # Discrete Event Simulation in C, pp 87-88.
  218.  
  219.     y = expovariate(alpha)
  220.     z = expovariate(1.0/beta)
  221.     return z/(y+z)
  222.  
  223. # -------------------- Pareto --------------------
  224.  
  225. def paretovariate(alpha):
  226.     # Jain, pg. 495
  227.  
  228.     u = random()
  229.     return 1.0 / pow(u, 1.0/alpha)
  230.  
  231. # -------------------- Weibull --------------------
  232.  
  233. def weibullvariate(alpha, beta):
  234.     # Jain, pg. 499; bug fix courtesy Bill Arms
  235.  
  236.     u = random()
  237.     return alpha * pow(-log(u), 1.0/beta)
  238.  
  239. # -------------------- test program --------------------
  240.  
  241. def test(N = 200):
  242.     print 'TWOPI         =', TWOPI
  243.     print 'LOG4          =', LOG4
  244.     print 'NV_MAGICCONST =', NV_MAGICCONST
  245.     print 'SG_MAGICCONST =', SG_MAGICCONST
  246.     test_generator(N, 'random()')
  247.     test_generator(N, 'normalvariate(0.0, 1.0)')
  248.     test_generator(N, 'lognormvariate(0.0, 1.0)')
  249.     test_generator(N, 'cunifvariate(0.0, 1.0)')
  250.     test_generator(N, 'expovariate(1.0)')
  251.     test_generator(N, 'vonmisesvariate(0.0, 1.0)')
  252.     test_generator(N, 'gammavariate(0.5, 1.0)')
  253.     test_generator(N, 'gammavariate(0.9, 1.0)')
  254.     test_generator(N, 'gammavariate(1.0, 1.0)')
  255.     test_generator(N, 'gammavariate(2.0, 1.0)')
  256.     test_generator(N, 'gammavariate(20.0, 1.0)')
  257.     test_generator(N, 'gammavariate(200.0, 1.0)')
  258.     test_generator(N, 'gauss(0.0, 1.0)')
  259.     test_generator(N, 'betavariate(3.0, 3.0)')
  260.     test_generator(N, 'paretovariate(1.0)')
  261.     test_generator(N, 'weibullvariate(1.0, 1.0)')
  262.  
  263. def test_generator(n, funccall):
  264.     import time
  265.     print n, 'times', funccall
  266.     code = compile(funccall, funccall, 'eval')
  267.     sum = 0.0
  268.     sqsum = 0.0
  269.     smallest = 1e10
  270.     largest = -1e10
  271.     t0 = time.time()
  272.     for i in range(n):
  273.         x = eval(code)
  274.         sum = sum + x
  275.         sqsum = sqsum + x*x
  276.         smallest = min(x, smallest)
  277.         largest = max(x, largest)
  278.     t1 = time.time()
  279.     print round(t1-t0, 3), 'sec,', 
  280.     avg = sum/n
  281.     stddev = sqrt(sqsum/n - avg*avg)
  282.     print 'avg %g, stddev %g, min %g, max %g' % \
  283.           (avg, stddev, smallest, largest)
  284.  
  285. if __name__ == '__main__':
  286.     test()
  287.